Differential peak-gene linkage analysis
Jack Collora
2023-06-29
Objectives
I want to determine if there is any change in the peak-gene linkage for differentially expressed genes across conditions. While the target genes are going to be primarily from the cytotoxic data, I do want to check all the populations to see if its a more common thing.
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(Signac))
results<-readRDS("~/gibbs/DOGMAMORPH/Ranalysis/Objects/20230628Linkpeakobj.rds")making the plotting function which will look at each of the individual conditions separately
#specifically made this for my morphine stuff, if someone else uses it in the future they should pay attention to the number of plots and the order. There are ways to make it more general but I opted to make it as close to what I need as possible
PlotLinked<-function(objlist, gene,idents.plot=c(""), up=10000, down=10000, expression.assay="RNA", assay="ATAC" ){
if(idents.plot[[1]]==""){idents.plot<-names(table(Idents(objlist[[1]])))}
p1<-CoveragePlot(objlist[[3]], region=gene, features=gene, extend.downstream = down, extend.upstream = up, assay = assay,expression.assay = expression.assay, idents = idents.plot)
p2<-CoveragePlot(objlist[[4]], region=gene, features=gene, extend.downstream = down, extend.upstream = up, assay = assay,expression.assay = expression.assay, idents = idents.plot)
p3<-CoveragePlot(objlist[[1]], region=gene, features=gene, extend.downstream = down, extend.upstream = up, assay = assay,expression.assay = expression.assay, idents = idents.plot)
p4<-CoveragePlot(objlist[[2]], region=gene, features=gene, extend.downstream = down, extend.upstream = up, assay = assay,expression.assay = expression.assay, idents = idents.plot)
p5<-CoveragePlot(objlist[[5]], region=gene, features=gene, extend.downstream = down, extend.upstream = up, assay = assay,expression.assay = expression.assay, idents = idents.plot)
p6<-CoveragePlot(objlist[[6]], region=gene, features=gene, extend.downstream = down, extend.upstream = up, assay = assay,expression.assay = expression.assay, idents = idents.plot)
return((p1+p2)/(p3+p4)/(p5+p6)+ patchwork::plot_layout(nrow = 3, byrow = FALSE) )
}
PlotLinked(results, gene = "CD4")## Warning: Removed 85 rows containing missing values (`geom_segment()`).
## Removed 85 rows containing missing values (`geom_segment()`).
## Removed 85 rows containing missing values (`geom_segment()`).
## Removed 85 rows containing missing values (`geom_segment()`).
## Removed 85 rows containing missing values (`geom_segment()`).
## Removed 85 rows containing missing values (`geom_segment()`).
PlotLinked(results,"OPRM1", up = 1e5, down = 1e5 )## Warning: Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
PlotLinked(results,"OGFR", up = 1e5, down = 1e5 )## Warning: Removed 13 rows containing missing values (`geom_segment()`).## Warning: Removed 13 rows containing missing values (`geom_segment()`).
## Removed 13 rows containing missing values (`geom_segment()`).
## Removed 13 rows containing missing values (`geom_segment()`).
## Removed 13 rows containing missing values (`geom_segment()`).
## Removed 13 rows containing missing values (`geom_segment()`).
selecting a few key genes that are upregulated paritcularly in the cyotoxic T and looking in just those populations
PlotLinked(results,"IFI44L", up = 1e5, down = 1e5, idents.plot = c("Effector_Th1_CD4_and_CD8","CD8_TEMRA") )## Warning: Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
PlotLinked(results,"PLCG2", up = 1e5, down = 1e5 , idents.plot = c("Effector_Th1_CD4_and_CD8","CD8_TEMRA") )## Warning: Removed 29 rows containing missing values (`geom_segment()`).## Warning: Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
PlotLinked(results,"BRD1", up = 1e5, down = 1e5 , idents.plot = c("Effector_Th1_CD4_and_CD8","CD8_TEMRA") )## Warning: Removed 5 rows containing missing values (`geom_segment()`).## Warning: Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
PlotLinked(results,"AKT3", up = 1e5, down = 1e5 , idents.plot = c("Effector_Th1_CD4_and_CD8","CD8_TEMRA") )## Warning: Removed 60 rows containing missing values (`geom_segment()`).## Warning: Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
PlotLinked(results,"CCL5", up = 1e5, down = 1e5 , idents.plot = c("Effector_Th1_CD4_and_CD8","CD8_TEMRA") )## Warning: Removed 22 rows containing missing values (`geom_segment()`).## Warning: Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
PlotLinked(results,"SPOCK2", up = 1e5, down = 1e5 , idents.plot = c("Effector_Th1_CD4_and_CD8","CD8_TEMRA") )## Warning: Removed 72 rows containing missing values (`geom_segment()`).## Warning: Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
PlotLinked(results,"NKG7", up = 1e5, down = 1e5 , idents.plot = c("Effector_Th1_CD4_and_CD8","CD8_TEMRA") )## Warning: Removed 3 rows containing missing values (`geom_segment()`).## Warning: Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
PlotLinked(results,"IFI44L", up = 1e5, down = 1e5)## Warning: Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
## Removed 11 rows containing missing values (`geom_segment()`).
PlotLinked(results,"PLCG2", up = 1e5, down = 1e5)## Warning: Removed 29 rows containing missing values (`geom_segment()`).## Warning: Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
## Removed 29 rows containing missing values (`geom_segment()`).
PlotLinked(results,"BRD1", up = 1e5, down = 1e5)## Warning: Removed 5 rows containing missing values (`geom_segment()`).## Warning: Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
## Removed 5 rows containing missing values (`geom_segment()`).
PlotLinked(results,"AKT3", up = 1e5, down = 1e5)## Warning: Removed 60 rows containing missing values (`geom_segment()`).## Warning: Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
## Removed 60 rows containing missing values (`geom_segment()`).
PlotLinked(results,"CCL5", up = 1e5, down = 1e5)## Warning: Removed 22 rows containing missing values (`geom_segment()`).## Warning: Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
## Removed 22 rows containing missing values (`geom_segment()`).
PlotLinked(results,"SPOCK2", up = 1e5, down = 1e5)## Warning: Removed 72 rows containing missing values (`geom_segment()`).## Warning: Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
## Removed 72 rows containing missing values (`geom_segment()`).
PlotLinked(results,"NKG7", up = 1e5, down = 1e5 )## Warning: Removed 3 rows containing missing values (`geom_segment()`).## Warning: Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
## Removed 3 rows containing missing values (`geom_segment()`).
devtools::session_info()## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1## - Session info ---------------------------------------------------------------
## setting value
## version R version 4.2.0 (2022-04-22)
## os Red Hat Enterprise Linux 8.8 (Ootpa)
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate C
## ctype C
## tz Etc/UTC
## date 2023-06-30
## pandoc 3.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/tools/ (via rmarkdown)
##
## - Packages -------------------------------------------------------------------
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [2] CRAN (R 4.2.0)
## BiocGenerics 0.44.0 2022-11-01 [1] Bioconductor
## BiocParallel 1.32.6 2023-03-17 [1] Bioconductor
## Biostrings 2.66.0 2022-11-01 [1] Bioconductor
## bitops 1.0-7 2021-04-24 [2] CRAN (R 4.2.0)
## bslib 0.4.2 2022-12-16 [1] CRAN (R 4.2.0)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.2.0)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.0)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.0)
## cluster 2.1.4 2022-08-22 [2] CRAN (R 4.2.0)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.2.0)
## colorspace 2.1-0 2023-01-23 [2] CRAN (R 4.2.0)
## cowplot 1.1.1 2020-12-30 [2] CRAN (R 4.2.0)
## crayon 1.5.2 2022-09-29 [2] CRAN (R 4.2.0)
## data.table 1.14.8 2023-02-17 [2] CRAN (R 4.2.0)
## DBI 1.1.3 2022-06-18 [2] CRAN (R 4.2.0)
## deldir 1.0-6 2021-10-23 [2] CRAN (R 4.2.0)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.2.0)
## digest 0.6.31 2022-12-11 [2] CRAN (R 4.2.0)
## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.2.0)
## ellipsis 0.3.2 2021-04-29 [2] CRAN (R 4.2.0)
## evaluate 0.20 2023-01-17 [2] CRAN (R 4.2.0)
## fansi 1.0.4 2023-01-22 [2] CRAN (R 4.2.0)
## farver 2.1.1 2022-07-06 [2] CRAN (R 4.2.0)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.0)
## fastmatch 1.1-3 2021-07-23 [2] CRAN (R 4.2.0)
## fitdistrplus 1.1-8 2022-03-10 [2] CRAN (R 4.2.0)
## fs 1.6.1 2023-02-06 [2] CRAN (R 4.2.0)
## future 1.32.0 2023-03-07 [1] CRAN (R 4.2.0)
## future.apply 1.10.0 2022-11-05 [1] CRAN (R 4.2.0)
## generics 0.1.3 2022-07-05 [2] CRAN (R 4.2.0)
## GenomeInfoDb 1.34.9 2023-02-02 [1] Bioconductor
## GenomeInfoDbData 1.2.9 2023-03-17 [1] Bioconductor
## GenomicRanges 1.50.2 2022-12-16 [1] Bioconductor
## ggforce 0.4.1 2022-10-04 [1] CRAN (R 4.2.0)
## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.2.0)
## ggrepel 0.9.3 2023-02-03 [1] CRAN (R 4.2.0)
## ggridges 0.5.4 2022-09-26 [1] CRAN (R 4.2.0)
## globals 0.16.2 2022-11-21 [1] CRAN (R 4.2.0)
## glue 1.6.2 2022-02-24 [2] CRAN (R 4.2.0)
## goftest 1.2-3 2021-10-07 [2] CRAN (R 4.2.0)
## gridExtra 2.3 2017-09-09 [2] CRAN (R 4.2.0)
## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.2.0)
## highr 0.10 2022-12-22 [1] CRAN (R 4.2.0)
## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.2.0)
## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.2.0)
## httpuv 1.6.9 2023-02-14 [1] CRAN (R 4.2.0)
## httr 1.4.5 2023-02-24 [1] CRAN (R 4.2.0)
## ica 1.0-3 2022-07-08 [2] CRAN (R 4.2.0)
## igraph 1.4.2 2023-04-07 [1] CRAN (R 4.2.0)
## IRanges 2.32.0 2022-11-01 [1] Bioconductor
## irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.2.0)
## jquerylib 0.1.4 2021-04-26 [2] CRAN (R 4.2.0)
## jsonlite 1.8.4 2022-12-06 [2] CRAN (R 4.2.0)
## KernSmooth 2.23-20 2021-05-03 [2] CRAN (R 4.2.0)
## knitr 1.42 2023-01-25 [1] CRAN (R 4.2.0)
## labeling 0.4.2 2020-10-20 [2] CRAN (R 4.2.0)
## later 1.3.0 2021-08-18 [2] CRAN (R 4.2.0)
## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.2.0)
## lazyeval 0.2.2 2019-03-15 [2] CRAN (R 4.2.0)
## leiden 0.4.3 2022-09-10 [1] CRAN (R 4.2.0)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.0)
## listenv 0.9.0 2022-12-16 [2] CRAN (R 4.2.0)
## lmtest 0.9-40 2022-03-21 [2] CRAN (R 4.2.0)
## magrittr 2.0.3 2022-03-30 [2] CRAN (R 4.2.0)
## MASS 7.3-59 2023-04-21 [1] CRAN (R 4.2.0)
## Matrix 1.5-4 2023-04-04 [1] CRAN (R 4.2.0)
## matrixStats 0.63.0 2022-11-18 [2] CRAN (R 4.2.0)
## memoise 2.0.1 2021-11-26 [2] CRAN (R 4.2.0)
## mime 0.12 2021-09-28 [2] CRAN (R 4.2.0)
## miniUI 0.1.1.1 2018-05-18 [2] CRAN (R 4.2.0)
## munsell 0.5.0 2018-06-12 [2] CRAN (R 4.2.0)
## nlme 3.1-162 2023-01-31 [1] CRAN (R 4.2.0)
## parallelly 1.35.0 2023-03-23 [1] CRAN (R 4.2.0)
## patchwork 1.1.2 2022-08-19 [1] CRAN (R 4.2.0)
## pbapply 1.7-0 2023-01-13 [1] CRAN (R 4.2.0)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.0)
## pkgbuild 1.4.0 2022-11-27 [1] CRAN (R 4.2.0)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.2.0)
## pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.2.0)
## plotly 4.10.1 2022-11-07 [1] CRAN (R 4.2.0)
## plyr 1.8.8 2022-11-11 [1] CRAN (R 4.2.0)
## png 0.1-8 2022-11-29 [1] CRAN (R 4.2.0)
## polyclip 1.10-4 2022-10-20 [1] CRAN (R 4.2.0)
## prettyunits 1.1.1 2020-01-24 [2] CRAN (R 4.2.0)
## processx 3.8.1 2023-04-18 [1] CRAN (R 4.2.0)
## profvis 0.3.8 2023-05-02 [1] CRAN (R 4.2.0)
## progressr 0.13.0 2023-01-10 [1] CRAN (R 4.2.0)
## promises 1.2.0.1 2021-02-11 [2] CRAN (R 4.2.0)
## ps 1.7.5 2023-04-18 [1] CRAN (R 4.2.0)
## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.2.0)
## R6 2.5.1 2021-08-19 [2] CRAN (R 4.2.0)
## RANN 2.6.1 2019-01-08 [2] CRAN (R 4.2.0)
## RColorBrewer 1.1-3 2022-04-03 [2] CRAN (R 4.2.0)
## Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.0)
## RcppAnnoy 0.0.20 2022-10-27 [1] CRAN (R 4.2.0)
## RcppRoll 0.3.0 2018-06-05 [2] CRAN (R 4.2.0)
## RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.2.0)
## remotes 2.4.2 2021-11-30 [2] CRAN (R 4.2.0)
## reshape2 1.4.4 2020-04-09 [2] CRAN (R 4.2.0)
## reticulate 1.28 2023-01-27 [1] CRAN (R 4.2.0)
## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.0)
## rmarkdown 2.22 2023-06-01 [1] CRAN (R 4.2.0)
## ROCR 1.0-11 2020-05-02 [2] CRAN (R 4.2.0)
## Rsamtools 2.14.0 2022-11-01 [1] Bioconductor
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.0)
## Rtsne 0.16 2022-04-17 [2] CRAN (R 4.2.0)
## S4Vectors 0.36.2 2023-02-26 [1] Bioconductor
## sass 0.4.5 2023-01-24 [1] CRAN (R 4.2.0)
## scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)
## scattermore 0.8 2022-02-14 [1] CRAN (R 4.2.0)
## sctransform 0.3.5 2022-09-21 [1] CRAN (R 4.2.0)
## sessioninfo 1.2.2 2021-12-06 [2] CRAN (R 4.2.0)
## Seurat * 4.3.0 2022-11-18 [1] CRAN (R 4.2.0)
## SeuratObject * 4.1.3 2022-11-07 [1] CRAN (R 4.2.0)
## shiny 1.7.4 2022-12-15 [1] CRAN (R 4.2.0)
## Signac * 1.9.0 2022-12-08 [1] CRAN (R 4.2.0)
## sp 1.6-0 2023-01-19 [1] CRAN (R 4.2.0)
## spatstat.data 3.0-1 2023-03-12 [1] CRAN (R 4.2.0)
## spatstat.explore 3.1-0 2023-03-14 [1] CRAN (R 4.2.0)
## spatstat.geom 3.1-0 2023-03-12 [1] CRAN (R 4.2.0)
## spatstat.random 3.1-4 2023-03-13 [1] CRAN (R 4.2.0)
## spatstat.sparse 3.0-1 2023-03-12 [1] CRAN (R 4.2.0)
## spatstat.utils 3.0-2 2023-03-11 [1] CRAN (R 4.2.0)
## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.2.0)
## stringr 1.5.0 2022-12-02 [1] CRAN (R 4.2.0)
## survival 3.5-5 2023-03-12 [1] CRAN (R 4.2.0)
## tensor 1.5 2012-05-05 [2] CRAN (R 4.2.0)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.2.0)
## tidyr 1.3.0 2023-01-24 [1] CRAN (R 4.2.0)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.0)
## tweenr 2.0.2 2022-09-06 [1] CRAN (R 4.2.0)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.0)
## usethis 2.1.6 2022-05-25 [1] CRAN (R 4.2.0)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.0)
## uwot 0.1.14 2022-08-22 [1] CRAN (R 4.2.0)
## vctrs 0.6.2 2023-04-19 [1] CRAN (R 4.2.0)
## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.2.0)
## withr 2.5.0 2022-03-03 [2] CRAN (R 4.2.0)
## xfun 0.39 2023-04-20 [1] CRAN (R 4.2.0)
## xtable 1.8-4 2019-04-21 [2] CRAN (R 4.2.0)
## XVector 0.38.0 2022-11-01 [1] Bioconductor
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.0)
## zlibbioc 1.44.0 2022-11-01 [1] Bioconductor
## zoo 1.8-12 2023-04-13 [1] CRAN (R 4.2.0)
##
## [1] /gpfs/gibbs/project/ya-chi_ho/jac369/R/4.2
## [2] /vast/palmer/apps/avx2/software/R/4.2.0-foss-2020b/lib64/R/library
##
## ------------------------------------------------------------------------------